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ABSTRACT 


A six degree of freedom robot manipulator arm, a PUMA 560, is calibrated using 
random subsets of available experimental calibration data. Some of these subsets produce 
good calibration results motivating the search for an optimum calibration procedure which 
will use a small number of poses. Statistical analysis of the joint excursions and end 
effector position variation in both "good" and "bad" subsets of poses were conducted. No 
significant statistical differences between them was discovered. The condition number of 
the Jacobian matrix is investigated as a potential measure of the accuracy which may be 
obtained from the subset under consideration. The condition number thus obtained 
contained too much variability to be a reliable predictor of accuracy. A computer 
simulation was conducted using a numerical optimizer to select the joint angles to be used 
for calibration. The optimizer studies failed to find an optimum set of poses for 
calibration. The conclusion of these studies is that there is no optimum set of poses to be 
used for calibration. An alternative hypothesis, that the resultant calibration accuracy 


depends only upon the accuracy of the measurements taken, seems to be proven. 
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I. INTRODUCTION 

The objective of this thesis was to investigate an optimum method of calibration of 
robotic manipulators. Calibration of manipulators seeks to improve their accuracy. 
Accuracy is measured in terms of both position and orientation (called the pose) of the 
manipulator end effector. It is the difference between the commanded pose and the 
achieved pose of the end effector. The achieved pose of the manipulator 1s a function of 
fixed geometric properties, such as link lengths, and of variable geometric properties, such 
as the angular displacement in revolute joints. A kinematic model is developed for the 
manipulator using both the fixed and variable geometric data. Errors between the pose 
predicted by the model and the pose measured in the laboratory for a typical manipulator 
have been determined experimentally to be 10 mm or more [Ref. 1]. These errors 
arise due to the differences between the nominal values and as built values of the 
geometric properties of the manipulator. Improving the accuracy requires a method of 
accurately determining these parameters. Several practical methods of calibrating 
manipulators have been investigated [Ref. 2]. 

Repeatability is another performance measure of a manipulator. Repeatability is the 
average measure of how closely the manipulator can achieve a pose which has been 
previously taught. Experiments have shown that the repeatability of a typical manipulator 


is on the order of 0.3 mm [Ref. 3]. Successful calibration of manipulators should 


be viewed as calibration which results in an accuracy which is close to the repeatability 
of the manipulator. 
Four basic steps in manipulator calibration have been identified and are briefly 
described as follows: 
¢ A closed chain kinematic model of the manipulator and measurement system is 
developed. During this process, identifiable parameters are determined and the 
measured quantity or quantities are specified. A set of error functions are derived 
from the difference in the measured quantities and the quantities predicted by the 
model. Nominal parameter values are provided by the manipulator manufacturing 
specifications, measurement system specifications and the location of the 
measurement system. 
¢ Next, experimental measurements are taken. These measurements are a function 
of the actual parameter values. Corresponding joint variable data is incorporated 
into the measurement set. 
¢ Identification of the parameters is performed utilizing the experimental data. This 
process consists of systematically adjusting the nominal parameters until the model 


predictions match the experimental data and hence the error functions become zero. 


¢ The final step involves incorporating the identified parameters into the software 
used to control the manipulator. [Ref. 4] 


In previous work, Swayze [Ref.5] performed the first three steps above on a 
PUMA 560 six degree of freedom manipulator arm and the fourth step was done in 
computer simulation. In his work, Swayze obtained calibration data consisting of 42 
poses and the joint angles associated with them. This data was arbitrarily separated into 
two groups of 21 poses. The first group was used as an experimental group and the 
second group was used as a control or reference group. Using a computer program, 
TEST, random sets of six poses were selected from the experimental group and a 


manipulator calibration was performed. The kinematic parameters identified by this 


calibration were used with the joint angles of the control group poses to calculate 
predicted poses using a forward kinematic solution. These predicted poses were then 
compared with the actual poses which had been experimentally obtained and an average 
position error was calculated. Swayze found that some of these randomly selected sets 
of six poses produced small position errors which were close to the repeatability of the 
PUMA 560 manipulator. The best set of six poses produced a position error of 0.46 mm. 
Figure (1) is a flow chart of how the program TEST operates. Swayze’s calibration data, 
which has not been previously published, is included as Appendix A. The program TEST 
was run for many iterations to investigate the range of position errors obtainable. The 
output files of program TEST were then screened by program SCREEN to select sets of 
poses which produced position errors less than 1.0 mm and greater than 20.0 mm. A 
summary of the output is tabulated in Appendix (B). These results raised several 


questions: 


¢ What is unique about the sets of poses which produced small position errors? 


¢ Can other small sets of poses be found which will also produce small position 
errors? 


¢ Is there an optimum way to select a small number of poses which will yield a 
calibration with accuracy approaching the repeatability of a manipulator? 


The answer to these questions is the topic of this thesis. 
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Figure 1. Flowchart for program TEST 


Il. THEORY 


A. KINEMATIC MODELING OF MANIPULATORS 
The theory presented and several diagrams of this chapter are based on material 
from Paul [Ref. 6]. The general arrangement of the material in this chapter 


closely follows the presentation of Swayze [Ref. 7]. 


1. General Coordinate System Transformations 

Robotic manipulators are constructed of multiple links connected by either 
revolute or prismatic joints. The kinematic model of the manipulator consists of a 
Cartesian Coordinate frame attached to each link with a set of transformation equations 
to describe positions in one coordinate frame in terms of another coordinate frame. 

Consider two coordinate frames, {0} and {1}, which have a common origin, 
where frame {1} is a produced from frame {0} by rotation of frame {0} by angle y about 
the x axis as shown in Figure (2). The position vector P may be represented in both 
frame {0} and frame {1} and in general will have different coordinates in each frame. 
Equations (1) are the transformation equations from frame {1} to frame {0} which can 
be readily seen by looking at a y-z planar projection of the two coordinate frames as 


shown in Figure (3). 





Y 
Xp : % 
°y, = ly,cosp - 'z,siny e 
P P » 
% =, Sind + 'z,cosy 
Expressing equations (1) in matrix format results in 
Oe 1 
*> 1 O 0 *p 
°y, = 10 cosy -siny 'Y (2) 
0 QO siny cos 1 
x, y y Za 
or °P =5R'P (3) 


where IR is the rotation matrix from frame {1} to frame {0}. Similar matrices can be 


derived for rotation about the y and z axes as well. For the remainder of this thesis, a 
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Figure 3. Coordinate frame rotation about the x axis projected onto the y-z plane. 


rotation about a coordinate axis will be represented as Rot(x,y), Rot(y,8), or Rot(z,6). 


The matrices that represent these rotations are given by equations (4-6). 


1 O 0 
(4) 


Rot(x,y) =|1 cosp -siny 
O sinW cosy 
cos8 QO sin®é 

Rot(y,8) =| 0 1 0 (5) 
-sin6 0 cos6 
cosh -sind 0 

Rot(z,d) =| sind cosh 0 (6) 
0 0 1 


Consider another coordinate frame, {2}, which is created by rotating frame 
{1} by an angle 0 about the y axis of frame {1}. In this case the vector P can be 


represented in frame {1} as 


'P =5R?P where 3R = Rot(Y,0) (7) 
and in frame {0} as 
°P =9R’P where °R = Rot(X,\) (3) 


or 
°P = "RRP =9R’P 
(9) 
where oR = OR sR. 
In the most general case, there can be translation as well as rotation of the 


coordinate frame. Let frame {1} be a translation from frame {0} by [x y z]' and frame 


{2} be a rotation of frame {1}. Then the position vector P will be expressed in frame 


{QO} as: 
x 
°p =!RP + ly (10) 
ré 


It would be convenient to have both translations and rotations from one coordinate frame 
to another represented as a single matrix transformation. To allow this requires use of 
a 4 x 4 matrix where the upper left 3 x 3 sub-matrix is the same rotation matrix as before, 
the mght hand column is the translation (x,y,z), and the last row of the matrix becomes 


[0 0 0 1]. This also requires augmented position vectors in the form 


(11) 


UD 


—_" 
= — 


where the 1 is a scale factor which will always be | in this work. Thus, the complete 


transformation of rotation and translation can now be expressed as: 


Nn 0. a. P, 


ip. |" Ov % ss (12) 
n. 0, a, P, 
1070 | 


where the n, 0, and a elements are the direction cosines of the x,y, and z axis of frame 
{1} (rotated frame) with respect to the original frame {0}. This 4 x 4 matrix is called 
the homogeneous transformation matrix and in general, represents a generalized rotation 
and translation of coordinate axes. If or is the transformation matrix which transforms 
position vectors in frame {1} to frame {0} then 
°p =iT'P. (13) 

If the position in the base frame {0} is known and the position in frame {1} is desired, 
then 

a ral n (14) 

= ol P. 

This presupposes that the inverse transformation matrix can be found. Rather than invert 


the T’ matrix, Paul [Ref. 8] describes the inverse as: 


Nn, ny n, —pn 


Tr =10 0, 0, —po (15) 
a, @, a, -pa 


where n, 0, a, and p are the column vectors of T and " - " represents the usual dot 


product of two vectors. 


Ze Roll, Pitch, Yaw and Euler Angle Transformations 

A given orientation may be specified in terms of three rotations. There are 
24 angle set conventions which may be used to specify these angles [Ref. 9]. In 
this work, rotations about fixed coordinate axes x-y-z will be used exclusively. This 
convention is referred to as roll, pitch and yaw (RPY) angles. The notation utilized will 
be 

RPY(,8,) = Rot(z,$)Rot(y,6)Rot(x,) Co 

which is interpreted as a rotation about x by an angle yw, followed by a rotation about y 
by an angle @ and lastly a rotation about z by an angle o. Miultiplying the rotation 


matrices for these 3 rotation together yields: 


cohc6 cohsOsy-shcy cohsOcy+sdsy 0 


RPY($,8,) = soc8 shsOsy+chcy sdhsOcy-cosy O (17) 
EEE -s8 cOsy cOcyy 0 
0 0 0 1 


where s and c have been used for sine and cosine for brevity. Multiplying equation (17) 
by the translation matrix, Trans(x,y,z), shown in equation (18), will yield a complete 


homogeneous transformation. Pre-multiplication would be used when the x, y and z 


10 


coordinates refer to the original coordinate axes prior to rotation and post-multiplication 


is used when the x, y and z coordinates are referenced to the rotated coordinate frame. 
10O0Ox 


0012 


0001 


Trans(x,y,Z) = 


Thus all transformations may thought of as a product of a rotation matrix and translation 
matrix of the form of equation (19). 

RPYT = RPY($,0,)Trans(x,y,z) (19) 
Figure (4) shows a generalized transformation with roll, pitch and yaw rotations with a 


translation of x, y and z in the new coordinate directions. 


3. Denavit-Hartenberg Transformations 

As previously stated, three rotations and three translations are required, in 
general, to transform one coordinate system to another. In robotic manipulators the 
attachment of successive links imposes constraints and therefore fewer parameters are 
required to fully describe the transformation. Although there are many different types of 
manipulators, they are all constructed of links attached at joints. Several systematic 
methods of attaching coordinate frames to links have been established. One such method 
is from Denavit and Hartenberg (DH) [Ref. 10] and is described below. 

Manipulator links are characterized by two parameters: link length a, and twist 
angle @, as shown in Figure (5). The link length a, is defined as the length of the 


common normal between the axes for joints n and n+l. The twist angle q, of a link is 
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Figure 4. Generalized Roll, Pitch, Yaw, Translation Transformation 


defined as the angle which joint axis n must be rotated about the common normal to 


reach joint axis n+l. A positive rotation is defined using the right hand rule considering 


the direction of travel from axis n to n+l to be the positive direction. Thus it can be seen 


that the link length must always be greater than or equal to zero while the twist angle has 


no such restriction. 
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The two links attached at a joint both have common normals as 


Joint n Joint n+l 





Angle 


On 
Link Length 


Figure 5. Link Length and Twist Angle 


described above. The distance between these common normals, measured along joint axis 
n, is defined to be the joint offset d,. The last parameter necessary to fully define the DH 
transformation is the joint angle 8,. The joint angle is defined to be the angle between 
the two common normals of a joint measured in a plane normal to the joint axis. 
Coordinate frames are established for each link of the manipulator by the 
following method. Link 0, which ts the base of the manipulator, is attached to link 1, the 
first movable link, via joint 1. Link 2 is attached to link 1 via joint 2 and so on. The 
origin of coordinate frame {n}, associated with link n, is located at the point of 
intersection of the common normal joining joint axes n and n+1 and the axis of joint n+l. 
In the case of intersecting joint axes the origin is chosen to be the point of intersection. 


The Z axis of frame {n} is chosen to be the axis of joint n+l. The x axis is chosen to 


be collinear with the common normal with the positive direction for x being in the 
direction of travel from joint axis n to joint axis n+l. In the case of intersecting joint 
axes the x axis is chosen perpendicular to the plane of intersection of the two axes. With 
the x and z axes thus established the Y axis is dictated by the use of mght hand 
coordinate systems. Note that the choice of the positive direction for z is arbitrary as is 
the positive direction of the x axis when the joint axes intersect. If joint axes are parallel 
then the position of the origin is chosen to make the link offset zero for the next link 
whose origin is defined. Figure (6) shows the assignment of link coordinate frames for 
a multiple link manipulator. 
With coordinate frames established for each link the relationship between link 

n-1 and link n is established by two rotations and two translations as follows: 

¢ rotate about z,, by an angle 90, which establishes the correct direction for x,; 

¢ translate along z,., a distance d, which gets the x, axis in the correct location; 


e translate along x, a distance a, which establishes the origin of frame {n} at the 
correct position; 


* rotate about x, by angle a, which establishes the correct orientation for z, and y,. 
These four steps are equivalent to equation (20) which is the DH transformation from 
frame {n-1} to frame {n}. Equation (21) is the matrix form of the DH transformation. 

"T = Rot(z,0) Trans(z,d) Trans(x,a) Rot(x,«) (20) 
To describe a position vector from the end effector frame in the base frame is as easy as 
multiplying the position P by the T matrix for each link of the manipulator as shown in 


equation (22). 
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Figure 6. Assignment of link parameters 9, d, a, and « 
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Modified Denavit-Hartenberg Transformations 


For ideal calibrations of manipulators it is important that the kinematic model 


used be proportionate. That is, a small change in one kinematic property should result 


in only small changes in the other kinematic properties. In most cases the DH method 


results in a proportionate model. However, in the case of parallel or nearly parallel axes 


the method becomes disproportionate. This is readily apparent from Figure (7) which 


He) 


Figure 7. Disproportionate affects of small axis variations on manipulator calibration 


represents two joints with nearly parallel intersecting axes. Let u-u’ represent the joint 
n-1 axis and w-w’ represent the joint n axis. A small rotation of axis w-w’ will cause a 
large change in the point of intersection. In the case of non intersecting nearly parallel 


joint axes a small rotation in one of them will cause a large change in the length of the 
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common normal a, between them. This small rotation may even cause non-parallel axes 
to become parallel such that no unique common normal exists. 

The modification to the standard DH transformations proposed by Hayati and 
Mirmirani [Ref. 11] will result in a proportionate model for consecutive revolute 
joints. The modification involves adding another rotation to the standard DH 
transformation and setting d, = 0. The following transformation can be used to arrive at 
either the standard DH or the modified DH (MDH) transformation. 

* T = Rot(z,0) Trans(z,d) Trans(x,a) Rot(x,«) Rot(y,B) (23) 


Equation (24), which expresses equation (23) in matrix form, becomes the standard DH 


cOcB-sOsasB -s6ca cOsBh+sBsacB acO 


n _ |SOcB+cOsasB cOca sOsB-cOsacB as6 (24) 
ce -casB Sa. cacB d| 
0 0 0 1 


transformation by setting B=0 and becomes the MDH transformation by setting d=0. All 
future reference to transformations between links of manipulators will be of the form of 
equation (24). This facilitates a standard computer code which is used for both types of 


transformations. 


5. Kinematic Chains and World Coordinate Frames 
When a series of links is joined together to form a robotic manipulator the 
series of transformations from the end effector coordinate frame to the base coordinate 
frame can be thought of as a kinematic chain. In this chain each transformation is 
equivalent to one link of the manipulator. In general the base frame {0} is internal to the 


manipulator and therefore knowledge of the position of the tool with respect to the base 


by 


frame is not very useful. A more useful reference coordinate frame would be the 
manipulator work space frame or world coordinate frame. This frame can be used to 
conduct measurements to be used for calibration and/or programming of the manipulator. 
Use of the world coordinate frame, designated frame {w}, requires use of an additional 
transformation to transform the base frame of the manipulator to world frame coordinates. 
This transformation may be done in one of two ways: 


¢ Use RPYT(0,9,w,x,y,z) to transform frame {0} to frame {w}. This adds six 
parameters to the kinematic model. 


¢ Use a DH style transformation to transform frame {0} to frame {w}. This requires 
only four parameters. 


At first glance, the second method appears to require fewer parameters but this is not the 
case. The fallacy is in the fact that the frame {0} for each method is different. In the 
first case frame {0} will be coincident with frame {1} with the result that two of the 
parameters of the MDH transformation from frame{Q} to frame {1} will be identically 
zero. In the second case four parameters will be required for the MDH transformation 
from frame {0} to frame {1}. Therefore, in both cases, a total of eight parameters will 
be required to transform between the world frame and the frame {1}. The last 
transformation in the chain is the sanctern aan between the tool frame and frame {n-1} 
for an n link manipulator. This transformation will be the RPYT transformation. The 


total kinematic chain is shown in equation (25). 


i My bey Ce (25) 
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6. Application to the PUMA 560 Manipulator 
The PUMA 560 manipulator consists of six links connected with revolute 
joints. The kinematic model is constructed using MDH transformations between the links 
and a RPYT transformation between link 5 and the end effector on link 6. The total 


number of parameters required for completeness is given by: 


N = 4R + 2P +6 (26) 


where N is the required number of independent parameters, R is the number of revolute 
joints, and P is the number of prismatic joints. [Ref. 12] The PUMA 560, with 
six revolute joints, thus requires 30 parameters for a complete model. A complete model 
is one which allows: 

¢ the reference (world) coordinate frame to be arbitrarily selected; 

¢ the zero position of each joint angle to be arbitrarily selected and; 

¢ the tool coordinate frame to be arbitrarily attached to link six of the manipulator. 
The MDH transformation has a specific definition of the zero position for each joint 
angle. To allow each joint zero position to be arbitrary a term 00, is introduced for each 
joint. Each 9, is the angle as measured by the joint encoder and 68, is the encoder offset 
from the zero position established by the MDH procedure. Using this definition, 9, for 
each joint is measured and considered fixed during a calibration while 59, is the parameter 
which is to be determined. 59, is considered to be constant for each joint throughout the 
range of travel of the joint angle. Table (1) shows the nominal kinematic parameters for 


the 30 parameter model. The values in parenthesis are defined to be zero (not part of the 
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model) and the underlined values represent the placement of the reference coordinate 
frame and tool frame used during this work. Figure (8) shows the assignment of 


coordinate frames for the PUMA 560 manipulator. 


B. PARAMETERS IDENTIFICATION METHODOLOGY 


1. The Numerical Optimizer ZXSSQ 
The IMSL subroutine ZXSSQ uses a Levenberg-Marquadt algorithm for the 
solution of non-linear least squares problems. The statement of the problem solved is 
given by equation (27). Each f'(x) is a user supplied objective function. The optimizer 
functions by calculating the gradient of the objective function at the current value of the 
design variable x using a finite difference approximation. It then changes x in a way that 


will reduce the objective function. This process is repeated until the objective function 


Table 1. PUMA 560 NOMINAL KINEMATIC PARAMETER TABLE 
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Figure 8. PUMA 560 Coordinate Frame Assignments using the DH Method 


minimum has been found or the user specified maximum number of iterations has been 
reached. ZXSSQ then passes the calculated value of the solution vector x back to the 


calling program. 


a), 


2. The Identification Jacobian Matrix 
For each calibration experiment conducted the kinematic parameters 
identification algorithm calculates an identification Jacobian matrix which relates the 
calibration poses and the kinematic parameters by the equation (28), 

SP = JéK (27) 
where OP is the difference between the experimentally measured pose and the calculated 
pose and OK is the difference between the calculated and the nominal kinematic 
parameters. J is the identification Jacobian matrix and will be calculated by the 
identification algorithm. In order to find 6K we must find the inverse of J. In general 
J will not be a square matrix so the equation must be inverted by post-multiplying both 
sides by the pseudo-inverse yielding the following equation for OK. 

8K = JT) TSP. (28) 
3. General Calibration Scheme 

Given the previously introduced 30 parameter model for the PUMA 560, 
which is based on the nominal kinematic parameter values, the actual values of the 
parameters may be found by performing a calibration. The pose of the end effector is 
measured and the associated manipulator joint angles are recorded. The joint angles are 
changed and the new pose is measured and joint angles recorded as before. This process 
is repeated until a sufficient quantity of data is obtained. Each pose consists of both 
orientation (3 rotation angles) and position (3 translations) so a sufficient number of poses 
is 30/6=5. There is, of course, some noise in the measurement system so a number of 


poses greater than the sufficient number is used so that the effect of the measurement 


i, 


noise may be reduced. For each pose and associated joint angles a forward kinematic 
solution is calculated based on the nominal model. This is the calculated pose referred 


to in the previous section. 


4. Calibration Implementation (Program ID6) 

The FORTRAN program ID6 is used to calculate the PUMA 560 manipulator 
kinematic parameters. The program reads in the nominal values of the parameters from 
a file (INPUT.DAT) and the calibration data from another file (POSE.DAT). The 
program then calls the IMSL subroutine ZXSSQ which calls a subroutine PUMA_ARM. 
PUMA_ARM calculates the forward kinematic solution for each set of joint angles using 
the current best guess of the kinematic parameters which have been supplied by ZXSSQ 
in the calling statement. It then computes the difference between the calculated and 
measured poses. These differences are returned to ZXSSQ as the objective function 
values. ZXSSQ then changes the kinematic parameters and calls PUMA_ARM again. 
When the objective function is reduced to approximately zero (that is, the measured and 
calculated poses are nearly identical) ZXSSQ has identified the kinematic parameters 
which it then passes back to the calling program. ID6 then writes these values to an 


output file and terminates. Figure (9) is a block diagram of ID6. 


23 


Read INPUT.DAT 
Read Calibration Data 


Output Identified 


Kinematic Parameters 





Calculate Pose Error 


Write Pose Error 


End 


Figure 9. Flow Chart for Program [D6 
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I. JOINT EXCURSION STUDY 


A. BASIS 

In this section the effect of total joint excursion during the conduct of a PUMA 560 
manipulator calibration is studied. The null hypothesis of this study is that the greater 
the total joint excursion in each joint of the manipulator, the better the resulting 
calibration accuracy. Since calibration attempts to accurately identify the position and 
orientation of revolute joints, each joint must be rotated through a range that is large 
enough to ensure that the joint axis may be properly identified. Consider a joint rotated 
through only a small joint rotation. Any calibration point measured during this rotation 
will be nearly collinear. When measurement noise is considered a large variation in the 
identified axis of rotation is possible. When calibration data is obtained for a larger 
rotation of the joint the variation in the uncertainty of the position of the axis of rotation 
is much smaller. Figure (10) illustrates the effect of small versus large rotation on the 


identified axis of rotation. [Ref. 13] 


B. METHOD 

As previously discussed in Chapter I, calibration data from the PUMA 560 (see 
Appendix A) was analyzed to find the average position error using six poses selected 
from the data set to conduct the calibration. The resulting accuracy versus pose numbers 


selected were tabulated in Appendix B. As was discussed in Chapter I, these results 
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Figure 10. Effect of small vs. large joint revolutions on joint axis identification. 


suggest that certain poses give better calibration than others. The study of this section 
uses the standard deviation of joint angles as a measure of the total joint excursion during 
each calibration experiment to see if there is a correlation between joint excursion and the 
resulting accuracy of calibration. It has been suggested [Ref. 14] that the 


accuracy of calibration is also a function of observation strategy and that for a fixed 


26 


number of observations the best calibration 1s obtained when the observations are 
randomly distributed over the manipulator work space. Therefore, both joint space and 
work space end effector excursions were compared with the resulting position accuracy 
of the calibration. The correlation coefficient between position error and each of these 
variables was calculated to see if a model could be constructed which would allow the 
prediction of the calibration accuracy obtainable. The method of calculating the standard 
deviations and correlation coefficients 1s described below. The standard deviations of 
joint space and end effector excursion were calculated using program STATI]. STATI 
used the output file created by Swayze’s program SCREEN, which consists of columns 
of position error with the pose numbers of the six poses used to calculate it, as an input 
file. STATI! reads the position error and pose numbers from this file and then reads the 
joint angles and end effector positions from Swayze’s calibration data files. The standard 
deviation of the joint excursion of each joint and the end effector excursion in x, y, and 
z 1s calculated and written to an output file along with the position error and the pose 
numbers. The output file was then input into MATLAB and the correlation coefficients 
between the position error and these standard deviations were calculated using 


MATLAB’s CORRCOEF command. 


C. RESULTS 
Table 2 shows the correlation coefficients of the joint angle standard deviations for 


each joint and work space position standard deviation versus position accuracy. 


2 


Table 2. POSITION ERROR CORRELATION COEFFICIENTS. 


Joint 5 0.0696 


Z Variation 













D. CONCLUSIONS 

The extremely low correlation coefficients for joint angle excursion and work space 
position excursion versus. position accuracy shows that there is no correlation between 
how much joint angle excursion occurs during a calibration and the resulting accuracy of 
the calibration. The same conclusion can be drawn about work space position excursion. 
Thus an alternative hypothesis, that calibration accuracy is independent of joint space 
excursion and work space excursion seems to be proved. If optimum poses for calibration 


exist they must be found through other methods. 


28 


IV. CONDITION NUMBER STUDY 


A. BASIS 
In the manipulator calibration equation, 

OP = JéK (29) 
introduced in Chapter II, 6P and J are not known exactly but have uncertainty. 
Measurement noise causes the uncertainty in oP and encoder noise causes the uncertainty 
in J [Ref. 15]. Consequently, there will be uncertainty in the solution for OK. 


If J is a square matrix, then 


_ IKI og yp IOI (30) 
K+ sxp Ow 
and 
ISK]. |6P|| (31) 
K+ skp WI 
where 
KV) = WI (32) 


is the condition number of the identification Jacobian matrix. In general, the manipulator 
calibration equation will be an over determined system and J will not be a square matrix. 
In this case, the condition number will be given by: 

KJ) = ITD MIG DL. e) 
It can be seen from equations (31) and (32), that for a small deviation in either the 


identification Jacobian or the measured poses for calibration, coupled with a large 


eo 


condition number, will cause a large uncertainty in 5K. That is, the accuracy of 
identification of the manipulator kinematic parameters will be degraded. Therefore it is 
hypothesized that the condition number of the identification Jacobian could be used as a 
predictor of the accuracy of the resultant calibration. 

Driels and Pathre [Ref. 16] have shown that use of the condition number 
appears to be a good predictor of the expected relative resultant accuracy when comparing 
calibrations conducted with the same number of observations but with different 
observation strategies. They achieved the same results using the condition number to 
differentiate between calibrations using the same observation strategy but with different 
numbers of observations. In this study, the utility of the condition number, as a predictor 
of calibration accuracy when comparing calibration experiments using the same 
observation strategy and number of observations, was studied. That is, can the condition 
number of the identification Jacobian be used to predict the best poses to use for 


calibration? 


B. METHOD 

The computer program TEST used by Swayze to calculate the position error from 
calibration data was modified to also calculate the L1 condition number for each set of 
poses tested. The program was run for 80 different sets of poses. The resulting 
accuracies and condition numbers were then compared. The poses selected for this study 
were in four groups as follows: 


e The twenty sets of poses which gave the lowest position error. The position error 
ranged from 0.464 to 0.619 mm. 
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¢ Twenty sets of poses which gave a position error in the range of 0.988 to 1.000 
mm. 


¢ Twenty sets of poses which gave a position error in the range of 20.15 to 21.64 
mm. 


¢ The twenty sets of poses which resulted in the largest position error of any sets 
considered. The position error ranged from 48.52 to 198.94 mm. 


C. RESULTS 
All of the position error versus condition numbers, along with the pose numbers 


used in the calibration, are tabulated in Appendix (B) and are summarized in Table (3). 


Table 3. POSITION ERROR VS. CONDITION NUMBER 


Position Condition 
Group Error (mm) Number(x 10°) 


l 0.464 - 0.619 0.095 - 1.8 


0.988 - 1.000 0.12 - 1.2 
20.15 - 21.64 3.2 - 125.4 


48.52 - 198.94 14.9 - >1000 





D. CONCLUSIONS 

Although the general trend of the condition number versus position error is as 
expected, (i.e. higher position error equates to higher condition number) the variability 
in condition number within a group of data is large enough to prevent any accurate 


predictions of the calibration accuracy given only the condition number. 
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V. ADS SEARCH FOR OPTIMUM POSES 


A. BASIS 

In Chapter I the motivation for a search for optimum poses for calibration was 
presented. In Chapters III and IV the condition number and joint excursion studies failed 
to detect any special significance between small sets of poses that produced superior 
calibration results and the sets which achieved poor calibration results. Clearly, if 
optimum poses exist, and they appear to, they must be found through other means. The 
hypothesis for the study of this chapter is that a numerical optimizer should be able to 
find the best set of poses to calibrate the PUMA 560 manipulator if an adequate objective 
function can be found. Since the goal in using an optimizer is to improve calibration 
accuracy it would seem reasonable to use position error as an objective function to be 
minimized. Since the absolute position error can never be known in a real manipulator, 
computer simulation studies were conducted so that its absolute position and orientation 


error could be known. 


B. METHOD 

A computer program, OPT6A, was written to conduct simulated calibration 
experiments. Figure (11) is a logical flow chart that shows how OPT6A works.r At the 
core of the program is a numerical optimizer, ADS, which is given an initial guess of the 


best sets of joint angles for conducting a calibration. In this simulation the exact 


a2 


Main Program Objective Function 


Initial Guess of 


Simulated Pose 


Joint Angles Measurement 





Kinematic Parameters 


_ Initialize ADS 
Identification 


Optimizer 





FWD Kinematic Sol’n 
to Calculate Poses 


Call Objective Func. 


Calc. Pos. & Orr. Err 
OBJ = POS+100*ORR 


Output Results 


Figure 11. Flow Chart for Program OPT6A 
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kinematic parameters for the manipulator to be calibrated are known. If the optimizer can 
find a set of joint angles which allows the exact kinematic parameters to be identified by 
the parameter identification algorithm, then this would be an optimum calibration. In 
each iteration of the optimizer the program calculates an objective function by a series 
of steps. 


¢ Simulated Pose Measurement Step: The program POSE [Ref. 17] is used 
to calculate the poses (associated with the joint angles) which would have been 
measured in the laboratory if this were a real vice simulated experiment. The poses 
produced have random noise injected to model the measurement process. 


¢ Parameter Identification Step: The poses produced above are used by the parameter 
identification algorithm similar to program ID6. The nominal kinematic parameters, 
with which the algorithm starts, are different by a known amount of 10 mm. for 
lengths and 5 degrees for angles. This ensures that the algorithm has a non-trivial 
identification problem. 


¢ Comparison Step: The program POSE is used again with the kinematic parameters 
identified in the previous step to find the forward kinematic solution of the joint 
angles in the verification data set. This set contains a relatively large number of 
sets of joint angles and poses that were produced with the exact kinematics and no 
measurement noise. The difference between the verification data set exact poses 
and calculated poses is then computed. The RMS position and orientation errors 
are then calculated. 


¢ Objective Function Construction: The orientation error has typically been found to 
be about 100 times smaller than the position error when positions are expressed in 
millimeters and angles are expressed in degrees. Therefore, the objective function 

used is the position error plus 100 times the orientation error. 
ADS then calculates the gradient of the objective function numerically and changes the 
joint angles in a direction which will reduce the objective function. This process repeats 


until ADS can find no further improvement in the objective function. The program then 


terminates and writes the joint angles and identified kinematic parameters to a file. 
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C. RESULTS 

The program was first tested with no simulated measurement noise. ADS conducted 
72 objective function calls and converged to the initial guess joint angle set. ADS takes 
two function calls per design variable to compute the gradient. The program was using 
6 poses with 6 joint angles each for a total of 36 design variables. Thus ADS converged 
to a solution which identified the exact kinematic parameters in the minimum number of 
function calls to achieve convergence. This was repeated several times with different set 
of randomly selected starting joint angles. The result was always the same, convergence 
in 72 function calls without changing the starting joint angles. 

The program was next run with simulated measurement noise. In this case the 
program did not converge. Due to the randomly applied simulated measurement noise 
the parameter identification algorithm did not identify the exact kinematic parameters. 
This caused a position and orientation error during each objective function call. 
Therefore, the objective function gradient is now non-zero and ADS changes the joint 
angles. With the new joint angles random noise is again applied in the simulated 
measurement step. The gradient calculated with these joint angles is also non-zero. After 
every gradient calculation (each 72 objective function calls) ADS changed the joint angles 
but they always stayed very close to the starting joint angles. This result was seen even 
if the starting point was a set of joint angles identified in Swayze’s work as giving a poor 


calibration. 


os 


D. CONCLUSIONS 

In the case of zero measurement noise the identification algorithm converges to the 
exact kinematic parameters on the first calculation of the gradient and this makes the 
gradient zero regardless of any dependency that the joint angles may have on resulting 
calibration accuracy. In the case of randomly injected simulation noise the identification 
algorithm could not identify the exact kinematic parameters and thus, the gradient can 
never be zero. However, the gradient due to the noise is small but large enough to keep 
the optimizer from converging. Therefore, if there is a dependency of calibration 
accuracy on joint angles (and hence the poses), the optimizer will tend to move the 
selected joint angles toward the optimum pose. No such movement of joint angles was 
detected so the conclusion must be that if such a dependency exists, it is only a very 
weak function of joint angles which produces such a small gradient that even the small 


simulated measurement noise completely obscures its effect. 


E. RETEST OF BEST POSES 

With the negative results of three totally separate investigations, serious doubts arise 
about the validity of the starting premise of this work, that is, that there exists select small 
sets of poses which will give nearly the same calibration accuracy as a large set of 
randomly picked poses. The research conducted thus far has found no evidence that this 
premise is true. Why then, did Swayze find such sets? 

In an attempt to answer this question an attempt was made to duplicate his results. 


Using the same PUMA 560 manipulator and the same coordinate measuring machine 
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(CMM) that Swayze used, the six poses previously found to give the smallest position 
error of any six poses were measured again. These were the poses identified in his data 
as pose numbers 8, 11, 17, 18, 19 and 20. The resultant six poses were input as the 
experimental group in program TEST with the resulting position error of 4.48 mm when 
compared to Swayze’s control group of poses. This error is approximately 100 times 
greater than the error of 0.464 mm originally obtained for these poses. The conclusion 
to be drawn from this result is that measurement error, though small, has a very big 


impact on the positional accuracy of the resultant calibration. 


F. POSSIBLE ALTERNATE HYPOTHESIS 

The results obtained in the previous section leads to the conclusion that the null 
hypothesis of this work must be rejected. A possible alternate hypothesis based on the 
no noise case described above is now proposed as follows. The positional accuracy of 
a manipulator which is achieved via calibration is a function of only the measurement 
noise in the calibration data and is totally independent of the joint angles used. This 
could explain why Swayze was able to find a small number of subsets of his calibration 
data which were able to give very good positional accuracy. The poses in these small 
subsets would be the ones which, for whatever reason, have smaller measurement error 
than the average pose. The poses in the small subsets which produced poor positional 


accuracy would be the ones which have a relatively larger measurement error. 
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VI. ZXSSQ SEARCH FOR OPTIMUM POSES 


A. BASIS 

In the last chapter, the program OPT6A which uses the optimizer ADS was 
presented. In this chapter, a similar program, OPT6C, is presented. OPT6C was written 
to use the IMSL subroutine ZXSSQ as the core subroutine vice the public domain 
program ADS. This was done for two reasons. First, there is technical support available 
for ZXSSQ while the user of ADS is not well supported. The second reason for rewriting 
OPT6A to utilize ZXSSQ in place of ADS is that OPT6A did not converge in the case 
where simulated measurement noise was injected into the calibration data. If the second 
optimizer, which uses completely different source code, also failed to find an optimum 


pose it would add to the body of evidence against the null hypothesis of this thesis. 


B. METHOD 

Figure (12) is a flow chart for OPT6C. The program uses two applications of 
ZXSSQ. The first application of ZXSSQ was renamed ZXSSQI1 to prevent fatal 
interaction between the two applications during program execution. The first application 
of ZXSSQ selects the joint angles for which a simulated calibration experiment is run. 
The selection is based upon trying to minimize the objective function. The objective 
function is calculated using the same four steps as previously described in Chapter V for 


program OPT6A briefly described as follows: 
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Figure 12. Flow Chart for Program OPT6C 
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¢ Simulated Pose Measurement 


¢ Parameter Identification: This step uses the other application of ZXSSQ in an ID6 
type program 


¢ Verification Step: Calculates position and orientation RMS errors 

¢ Calculate the objective function and pass it to ZXSSQI1. 
Based on the results of the simulated calibration experiment ZXSSQI changed the values 
of the joint angles and ran another simulated calibration experiment. When the objection 
for function for ZXSSQI1 does not change the program has converged to an optimum set 


of joint angles for calibration. 


C. RESULTS 

The results from OPT6C were exactly the same as for OPT6A for both the no 
measurement noise and measurement noise added cases. In the no simulated 
measurement noise case the outer loop optimizer, ZXSSQ1, converged after enough 
function calls to calculate the gradient to the starting joint angle set. When simulated 
measurement noise is added the program does not converge but does change the joint 
angles after each gradient evaluation is complete. If a gradient due to a dependence on 
joint angles were present the joint angles selected would tend to change in the direction 
of the optimum set of joint angles. Just as in program OPT6A no trend in the movement 
of the joint angles is observable and the joint angle set stays essentially the same as the 


Starting set. 
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D. CONCLUSIONS 

In the no noise case the optimizer converges to the starting joint angle set 
independent of the starting joint angles independent of which optimizer is used to select 
the joint angles. In the case where simulated measurement noise was added to the 
calibration data the noise affected the gradient computed at each step only a small amount 
but even this small amount totally obscured any dependence of positional accuracy on 
joint angles used. That is, the gradient due to the measurement noise is superimposed on 


an essentially flat topology. 
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Vil. RANDOM SEARCH FOR OPTIMUM POSES 


A. BASIS 

In the previous two chapters two numerical optimizers were shown to have been 
unable to detect a dependency of the objective function on joint angles. Two possible 
reasons for this are apparent. First, there may be no dependency to be found which 
means that there are no optimum poses to be found, or second, the dependency may be 
a very weak function of joint angles which the optimizers failed to detect because they 
have tested a neighborhood around the starting point which was too small. The purpose 
of the investigation of this section was to check large neighborhoods for objective 
function gradients. If an optimum set of poses exists it will result in an objective 


function at that point which has a lower value than other point. 


B. METHOD 

In this chapter the method of search for an optimum pose was changed from that 
of the previous two chapters. Instead of changing the joint angles based on the objective 
function, as was done in both programs OPT6A and OPT6C, in this search random points 
were picked in a neighborhood around the starting point. At each point selected the 
objective function is calculated. Each time the objective function is calculated the 
simulated measurement noise injected in the calibration poses is different because it is 


randomly selected. To get a measure of the range of calibration error resulting from this 
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noise the objective function is calculated 100 times at each randomly selected point in the 
neighborhood of the starting point. The mean and standard deviation of these 100 
objective function calculations 1s computed and written to a file. This process is repeated 
for different sized neighborhoods ranging from £5 to £180 degrees of variation from the 
starting angle for each of the six joints of the PUMA 560 manipulator. The objective 
function used is the same as the objective function used in OPT6A and OPT6C, position 
error plus 100 times orientation error. All of the data obtained for this chapter was 
obtained using a version of the program OPT6C that was modified slightly to defeat the 
joint angle selection subroutine ZXSSQI and call the objective function directly vice let 
it be called from ZXSSQI1. Each of the neighborhoods tested was tested with 100 random 
points. The 180 degree variability case was run ten times to give a total of 1000 random 


points tested throughout the working volume of the manipulator. 


C. RESULTS 

Each of the 20 random search data runs conducted used the same set of six poses 
and associated joint angles as the starting point (center of the neighborhood). The test 
results are presented in Table (4). The objective function of the starting point was 
calculated 100 times to give a comparison for each of the other points tested. The 
degrees column represents the magnitude of the maximum variation in each joint angle 
during the test. Each row in the table (except the 0 degree row which is neighborhood 
consisting of only one point) represents the low and high of 100 data points in the 


specified neighborhood of the starting point. The low and high refer to the 


43 


Table 4. OBJECTIVE FUNCTION DATA SUMMARY 


Low Data High Data 
N 
l 


= 


NA 
_ sor | 099 | 62s 
2.79 


186.41* 5.38* 
[| ae [as || 
fr 
130-8 | 484 | 105 197 


* The 186.41 mean was the highest with 13.28 next highest. It appears that this 
point was very near a singular point. 


5 
10 
20 
30 
30 


lowest/highest mean of the 100 points tested for that particular neighborhood. In each of 
the runs presented, the measurement error injected was 0.2 mm for position and 0.2 


degrees for angles times a random number between -1.0 and 1.0. 


D. CONCLUSIONS 

A review of the means and standard deviations of the low data shows that the mean 
objective function does not change appreciably throughout the working volume of the 
manipulator. This means that there is no optimum set of six poses (and associated six 
sets of joint angles) which will give a calibration accuracy that is as good as a large set 
of calibration data. That is, if the calibration data all contain a random measurement error 
there can be no guarantee of a good calibration unless a large number of poses is used 
to calibrate the manipulator. A review of the high data shows that there are, however, 
sets of poses which produce relative poor calibration results. These sets are those which 


in some way approach a singularity of the manipulator. 
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VIN. DISCUSSION, CONCLUSIONS AND RECOMMENDATIONS 


The null hypothesis of this work was that there exists select small subsets of poses 


which will give calibration accuracy approaching the repeatability of the manipulator. 


Previous experimental work conducted by Swayze had shown that there were such small 


subsets. The goal of this thesis was to investigate these subsets and discover what is 


unique about them, how to find them and then use them to predict calibration accuracy 


prior to actually conducting the calibration. In pursuit of this goal a three pronged 


investigation of the problem was conducted as follows: 


Investigate the condition number of the calibration Jacobian matrix. 


Conduct statistical analysis of sets of poses which yielded good and poor calibration 
results. 


Conduct computer simulation studies searching for optimum poses. 


During these investigations several interesting discoveries were made: 


The condition number of the calibration Jacobian matrix is not an effective 
predictor of calibration accuracy. 


There is no Statistical difference between either the joint space excursions or work 
space position excursions of subsets of calibrations which produced good calibration 
accuracy and those which produced poor calibration accuracy. 


The computer simulation studies could not find a set of poses which produced a 
better calibration than a similar number of poses chosen at random. 


One set of poses (obtained by Swayze) which were known to give a calibration 


accuracy approaching the repeatability of the PUMA 560 manipulator gave a much 
poorer result when they were measured again. 
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Based on these discoveries it is concluded that optimum poses for conducting kinematic 
parameters identification do not exist. Thus, the following alternate hypothesis appears 
to be proven. The accuracy which results from conducting a manipulator calibration is 
a function of the pose measurement error only. It has been shown that without 
measurement noise any set of poses used during calibration wul yield the exact kinematic 
parameters. 
A recommendation for further study is to investigate an optimum way to overcome 

the affects of measurement error. Given that there will be measurement noise: 

¢ Is the measurement noise randomly distributed? 

¢ Does the measurement noise also have a bias? 


¢ (Can a method be discovered to tell which poses have the most noise? 
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APPENDIX A 
The following 42 data sets each consist of six joint angles and a four by four 
transformation matrix. The last row of the transformation matrix is [0 O O 1] and is 
inserted by the computer program which analyzed the data, TEST.FOR. This data was 
obtained by Swayze [Ref. 18] and is published here only because it has not 
previously been published and the data has been used in the research for this thesis. The 


numbers in parenthesis are the pose numbers used Swayze. 


Experimental Group of Poses from file 121.dat: 


(1) 

108.6970 39..33100 ayvelZ 200 169.9640 ~99 .93700 155.5300 
O..77 143 O.6178214 ~0.1544949 752.0477 
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0.5205824 0..8528200 2-1506701E-02 330.3344 
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-109 .6930 137-0600 -170.7660 3.812009 =90 .80800 265.9100 
-“V-tolesos =0:.:9553381 =0.2543225 Zo. 16502 

0.818687 I 2.1206835E-02 -0.5741130 433.8619 

05525395 ~Q3 2958570 Hed 792214 205.5811 

(4) 
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0.3266630 0.4924835 0.8074185 37 215068 
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(5) 

=35 60600 137.0000 153.6440 99 -15800 99.95400 12 32500 
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0.8433744 -0 .4885970 =022123972 200.5175 

-0.2607758 -0.7246256 0.6346267 279.1954 

(6) 
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Control Group of Poses from file 221.dat 
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APPENDIX B 
The following data is from two sources. The position error was generated by 
Swayze using the program TEST. The six pose numbers (which represent the poses from 
Appendix A) yielded a calibration with position errors shown. The L1 condition numbers 
were generated using a modified version of TEST which was rerun using the same 


calibration data. 


Group 1 Data: 


Position Lil Condition 

Eeror Pose Numbers Used Number 
0.464226 8 iB i: 18 19 20 489388. 
0.491716 3 8 9 Tek 18 19 IL ALS AAS uilees 
O-526177 Z 9 Ab al 7 18 Z0 1103019. 
O2S512113 Z 9 itil 13 17 20 610886. 
0.550705 3 6 8 iat 7 18 94999. 
0.577800 4 6 ill iy 18 19 312798. 
02580232 5 9 1! 17 18 20 1798603. 
0.581695 3 & 8 10 iL 18 207s 
0.5835 4 Z 3 5 qi) 9 lial 542206. 
0.584908 iL 9 13 18 20 Zl 10412637 
0.584987 2 5 9 Uae 18 20 Bias ee 
0.585944 2 9 acca 2 Ly 18 5579372 
0.606438 9 13 15 18 20 2611 582533. 
0.608049 2 4 4) 9 da 17 265467. 
0.608943 9 14 ig 18 19 20 790714. 
0.614340 3 6 8 7 18 19 961907. 
0.615247 4 a 10 16 iy 20 902923. 
O2622127 3 6 fi 9 17 19 235194. 
0.622639 5 9 Ls 1% 18 20 245892. 
0.619353 I 8 9 IS 7 20 744a2755 

Group 2 Data: 

Position Ll eeonditicn 

Error Pose Numbers Used Number 
0.988286 5 6 9 ee, 13 19 395734. 
0.988995 1 3 4 it ike 19 349937. 
0.989457 2 3 4 10 lea 17 10427 7 
0.989890 q 9 eZ ILS 19 Zk 32027 0: 
02991495 il 3 6 8 iz Ly PAPAS 0) 1b. 
0.992349 5 8 ital 14 183 19 486052. 
0.992768 6 ILA LS 16 18 19 709844. 
0.993289 Z 6 10 Le 19 Zi 289391. 
0.993329 a 8 9 cs 7 19 1025013. 
0.993746 7 ee 14 di 19 20 124090. 
0.993841 6 9 7 18 19 real 1S 26927 
0.994538 2 4 9 10 7 Ze 578332. 
0.995158 il 4 Z 9 ea eZ 716899. 
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Group 2 Data continued: 


Position 
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Group 3 Data: 
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Group 4 Data: 
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Number 
S919 12. 
380416. 
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